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“Hypothetical” heavy particles dynamics in LES 
of turbulent dispersed two-phase channel flow 

By M. Gorokhovski and A. Chtab f 


1. Motivation and objectives 

The extensive experimental study of dispersed two-phase turbulent flow in a vertical 
channel has been performed in Eaton’s research group in the Mechanical Engineering 
Department at Stanford University (see Kulick et al. (1994) and Fessler et al. (1994)). 
In Wang & Squires (1996), this study motivated the validation of LES approach with 
Lagrangian tracking of round particles governed by drag forces. While the computed 
velocity of the flow have been predicted relatively well, the computed particle velocity 
differed strongly from the measured one. Using Monte Carlo simulation of inter-particle 
collisions, the computation of Yamamoto et al. (2001) was specifically performed to model 
Eaton’s experiment. The results of Yamamoto et al. (2001) improved the particle velocity 
distribution. At the same time, Vance & Squires (2002) mentioned that the stochastic 
simualtion of inter-particle collisions is too expensive, requiring significantly more CPU 
resources than one needs for the gas flow computation. Therefore, the need comes to 
account for the inter-particle collisions in a simpler and still effective way. To present 
such a model in the framework of LES/Lagrangian particle approach, and to compare 
the calculated results with Eaton’s measurement and modeling of Yamamoto is the main 
objective of the present paper. 


2. Equation for the particle motion along smoothed trajectory 

In the high-Reynolds number flows, three random forces act on the motion of solid 
particles: (i) particle-small scale turbulence interaction (affecting mostly light particle 
dynamics); (ii) inter-particle collisions (important for relatively heavy particle motion; 
this randomness motivated the present work); (iii) particle- wall interaction (due to the 
wall roughness). Accounting for all details of this randomness in LES is formidable. The 
problem is how to describe effectively the random Lagrangian motion of solid particles. 
One of the possibilities is as follows: By analogy with kinetic theory, we may introduce 
a ‘hypothetical’ particle, which is moving along the smoothed trajectory averaged over 
random collision trajectories. The smoothed acceleration of such solid particle changes 
the local momentum in the gas phase, thereby affecting the acceleration of neighboring 
solid particles. To derive the equation of such a smoothed motion, let us consider the 
dispersed two-phase turbulent flow as a time evolution of a system of interacting stochas- 
tic fluid (Pope 2000) and solid particles. The kinetic description is then specified by the 
distribution function of particles with position 1?(xi, X 2 , £ 3 ) and velocity W(vi, V2, V3), 
/(T\ V,f), in the six-dimensional space, for solid particles, /^(Yp, "V p , t) and fluid par- 
ticles, fg(~r g ,^ g ,t). For solid particles, the Boltzman equation can be written in the 
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following form: 

= I(f p , f p ) + I(f p , f g ) (2.1) 

where I(f p ,f p ) and I(f p ,f g ) are operators for solid-solid and solid-liquid particle colli- 
sions, respectively. The first collision operator in (2.1) is usually modeled assuming pair 
hard sphere collisions with prescribed efficiency (see Vance & Squires (2002), Lavelieville 
et al. (1997), Simonin et al. (2002), Sommerfeld (2001), Wang et al. (2000), Mei & Hu 
(1999), Sigurgeirsson et al. (2001)). The form of the second operator is principally un- 
known. The total number of solid particles does not change due to any collisions, as well 
as the total momentum of solid particles is unchanged due to solid-solid collisions, i.e.: 

j d\ p I(f p J p ) = 0 (2.2) 

J d 3 x p I(f p J g ) = 0 (2.3) 

J d 3 v p I(f p , fp)v p = 0 (2.4) 


However 

J d 3 x p I(f p Jg)v p ± 0 


(2.5) 


Introducing the solid particles density, n p (in usual terms of kinetic approach) and the 
particle velocity averaged by f p , (v p )/ = d- / d 3 v p f p v p , let us assume that the solid- 
solid particle collision occurs frequently such that the correlation of fluctuations of solid 
particle velocities can be identified with the particle temperature of ordinary statistical 
mechanics, T p : 


( 


v — v k/3>/ p = o ( 


^P,Ct ' 


V-) fM s a 0 


( 2 . 6 ) 


Here /m is Maxwellian distribution, m p is the mass of the considered solid particle 
and 6 a /3 is the Kronecker delta. Multiplying (2.1) by the velocity of particle, integrating 
over all these velocities and using (2.6), one yields: 


rf ( v P,/3)/ P 

dt 



(2.7) 


where the first term implies the change of the ‘particle pressure’ along the smoothed 
trajectory of the particle and the second term is the mean rate of the particle velocity 
change due to particle-turbulence interaction. These two terms have to be modeled. 
The first term is modeled as follows. We introduce the ‘temperature’ of stochastic fluid 
particles, T tur g . According to its kinetic definition, one writes: 

T t ur, g = 2 - Pg VOL L (V*)/2 (2.8) 

where p g is the gas density, V OL l may be viewed as a volume over integral spatial scale 
and (v'g)/2 is the kinetic energy of turbulence in gas flow. Assuming a relaxation of the 
statistical temperature of solid particle to the statistical temperature of liquid one, one 
writes: 

= mur,g - T p ) (2.9) 

where /3 is an exchange frequency parameter. If (3 is very large, the statistical temperature 



‘‘ Hypothetical ” heavy particles dynamics 207 

of particle may be assumed to be locally in a ‘thermodynamic’ equilibrium with the 
surrounding turbulence: 

Tp = T t ur,g (2-10) 

Concerning the second term in (2.7), for heavy particles, it can be presented by the draft 
force: 

(( v S,/3)/ s - ( v P,/3>/p) / T P (2- 11 ) 

where t p is the Stokes time. With the particle Reynolds dependency from Clift et al. 
(1978), it writes: 

_ Pv d l 1 

p 18 p g v g 1 + 0.15 Re° p 6s7 1 ' ’ 

In this formulation, the equation (2.7) without the first term represents usual LES/ 
particle tracking procedure. 


3. Computation procedure 

LES of particle-laden flow is performed for conditions corresponding to the Eaton’s 
experiment. In this experiment, a set of measurements has been performed on a ver- 
tical fully-developed air-channel flow laden with spherical particles of different mass 
loading. From this measurement, we chose the case with particles of copper (density 
p p = 8800 kg /in 3 , diameter d p = 70 pm and 20% of mass loading). The numerical code 
developed at CTR, Stanford University by Pierce & Moin (2001), was adapted in this 
paper for IBM PC. In this code, the sub-grid momentum transport term was modeled by 
dynamic approach of Germano et al. (1990). In the present work, this code has been cou- 
pled in ‘two-ways’ with Lagrangian particle solver, according to (2.7)-(2.11). The second 
order Runge-Kutta scheme has been used for computation of the particle motion. 

To define the local statistical temperature of turbulence in (2.8), the volume was as- 
sociated with the control volume of the finite-difference mesh. Two different expressions 
for the ‘temperature’ of turbulence have been used to calculate (2.8): 

T t urb,g = \p g A x A y A z 5> s , q - v g , Q ) 2 / 2 (3.1) 

a 

and for the mean one 

(' Tturb.g ) = ^ PgA x A v A z ^ ((v 9iQ - v g>ct ) 2 ) /2 (3.2) 

ot 

The numerical algorithm of two-way momentum coupling was implemented. The par- 
ticle equations have been solved using the second order Runge-Kutta method. For the 
gas velocity at a particle position, the linear interpolation scheme has been used. The 
implementation of higher order schemes of interpolation did not give an explicit advan- 
tage, requiring at the same time a substantial computational effort. The computations 
were performed at Reynolds number based on friction velocity and channel half-width 
of 644 (corresponding to Reynolds numbers of 13800 based on centerline velocity and 
channel half- width). Parameters of the computation have been chosen similar to Wang 
& Squires (1996) with 64 x 65 x 64 grid points for the flow resolution in the x, y and z 
directions, respectively, that covered the computational domain hTiS/2 x 26 x tt5/2. In 
the streamwise and spanwise directions, uniform grid was used. In the direction normal 
to the wall, the non-uniform stretched grid has been used with first velocity position at 
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■y + = 0.88. Periodic boundary conditions were used for gas phase and for particles at free 
boundaries. 


4. Results 

For unladen and laden turbulent flow with different mass loading, Fig. la shows the 
mean stream-wise gas-velocity obtained in Eaton’s experiment and in the present com- 
putation. It has been noted in the experiment that particles do not change practically 
the mean velocity profile of gas. As it is seen from Fig. la, the computed mean gas ve- 
locity field is also not influenced by the presence of particles. Fig. lb and lc show the 
turbulence intensity in stream- wise and wall-normal directions. The Eaton’s experiment 
predicted an attenuation of the turbulence by particles. This is also seen from compu- 
tation in Fig. lb and lc. To demonstrate the strong effect of the turbulence attenuation, 
the centreline experimental and computed values of kinetic energy and viscous dissipa- 
tion are presented in Fig. 2. versus different mass loading. In Fig. 3, the profiles of mean 
particle velocity and of its variance are shown. Fig. 3a gives the particle mean streamwise 
velocity. It is clearly seen that similar to measurements and computations of Yamamoto 
et al. (2001), the present computation gives the flattened profile. In Fig. 3b, the computed 
r.m.s. streamwise particle velocity is compared with the measurements. The computation 
is in agreement with experiment except near the channel centre. At the same time, it is 
seen from Fig.3c, that computed r.m.s. of wall-normal velocity is overestimated against 
the measured values (the data from Yamamoto et al. (2001) show an underestimation of 
measurements) . 

5. Conclusion and future work 

The recent LES computation of particle-laden turbulent flow showed that the particle 
velocity distribution can be predicted relatively well, if the inter-particle collisions are 
included in the simulation. However the Monte Carlo simulation of inter-particle colli- 
sions, which is used in those computations, is too expensive for practical applications. In 
this paper, a simplified approach of two-phase flow simulation is proposed. By analogy 
with kinetic theory, we introduced a ‘hypothetical’ particle, which is moving along the 
smoothed trajectory averaged over random collision trajectories. An equation of such a 
smoothed motion has been derived, using hypothesis of a ‘thermodynamic’ equilibrium 
between the statistical temperature of the particle and the surrounding turbulence. The 
Lagrangian tracking of ‘hypothetical’ particles was performed along LES computation of 
a vertical fully-developed air-channel flow (experiment in Eaton’s research group). It has 
been shown that our computation is in agreement with Eaton’s experiment and compu- 
tation, where the inter-particle interaction has been simulated by hard-sphere collisions 
with prescribed efficiency. At the same time, to account for inter-particle interactions, 
the presented model does not require additional CPU time. Our future work concerns 
the further development of this approach for interaction between particles & rough wall. 
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Figure 1. Comparison of computed gas velocity profiles with measurements of Kulick et al. 
(1994) at different mass loading: a) streamwise mean velocity; b) r.m.s. of streamwise velocity; 
c) r.m.s. of wall-normal velocity. 
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Figure 2. Attenuation of gas-phase turbulent kinetic energy and dissipation with particle 
mass loading: filled symboles - computation; empty symbols - measurement. 
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Figure 3. Comparison of computed particle velocity profiles with measurements of Kulick et 
al. (1994) and results from Yamamoto et al. (2001) with inter-particle collisions: a) streamwise 
mean velocity; b) r.m.s. of streamwise velocity; c) r.m.s. of wall-normal velocity. 
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